Branching Processes and Evolution at the Ends of a Food Chain 



G. Caldarelli, C. Tebaldi 
S.I.S.S.A./I.S.A.S., V. Beirut 2-4, 1-34013 Grignano di Trieste (TS), Italy 
I.N.F.M., Istituto Nazionale di Fisica della Materia, Sezione di Trieste 

A. L. Stella 

INFM-Dipartimento di Fisica e Sezione INFN, Universita di Padova, 

1-35131 Padova, Italy 

Abstract 

In a critically self-organized model of punctuated equilibrium, boundaries 

determine peculiar scaling of the size distribution of evolutionary avalanches. 

This is derived by an inhomogeneous generalization of standard branching 

processes, extending previous mean field descriptions and yielding u = 1/2 

together with t' = 7/4, as distribution exponent of avalanches starting from 

species at the ends of a food chain. For the nearest neighbor chain one obtains 

numerically r' = 1.25 ±0.01, and Tf irst = 1.35 ±0.01 for the first return times 

of activity, again distinct from bulk exponents. 

PACS: 02.50.-r, 87.10.±e, 05.40.±j, S.I.S.S.A. Ref. 21/96/CM 
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Branching processes (BP) occur in many fields of physics and biology, ranging from 
nuclear reactors, to polymers and population dynamics Within the context of self- 
organized criticality (SOC), introduced by Bak, Tang and Wiesenfeld [[J, BP, or correlated 
versions of them, are expected to underlie the physics of many models, describing sandpiles 
earthquakes ]5|], river networks || or species mutations JTJ. By evolving long enough, 
these models self organize in stationary critical states with long-range correlations in space 
and time, and with avalanches of activity occurring at all scales. Avalanches are often 
believed to be described in terms of critical BP in the mean field (MF) limit. In the 
present Letter we introduce and solve an inhomogeneus generalization of the standard BP. 
This allows us to determine peculiar scaling properties of BP at boundaries. In a unifying 
perspective, such properties provide a substantial extension of previous MF descriptions of 
SOC models. 

Bak and Sneppen (BS) |7],§] introduced a SOC model describing an ecology of interacting 
species evolving by mutation and selection. This model provides an illustration of the 
mechanisms determining intermittency (punctuated equilibrium ||) and scaling |TI| in the 
evolutionary activity. Below we show that such intermittency and scaling have a richer 
structure than appreciated so far. Indeed, at the level of universal properties, it is possible 
to draw a clear cut distinction between evolutionary activities occurring in the "bulk" and 
at the "boundary" of an ecology. Bulk and boundary refer to different locations of a given 
species within the network of interactions with other species conditioning its evolution. 

In a coarse grained, simplified description, BS associate to the i-th species of an ecology 
a single fitness parameter, Xi, (0 < X{ < 1). Xj represents the ability of species i to survive: 
the higher x«, the higher the barrier to overcome in order to switch a mutation in the species. 
A genetic mutation changes the barrier of the species and modifies also the barriers of the 
other species interacting directly with it. This interaction should represent the fact that two 
species, e.g., take part in the same food chain. Sites of a lattice can be used to represent 
the species: in this case nearest neighbor (n.n.) species can be assumed as directly related 
biologically, and thus interacting. 
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The dynamical evolution rules are as follows. Starting from an initial fitness landscape, 
the i with lowest x, i m in, is selected to undergo a mutation and its fitness x imin is modified 
into a new one, chosen at random. Due to the interaction, also some neighbors of i m i n get 
modified x's, as effect of the previous mutation. For a linear chain with n.n. interactions 
this implies that x,- _i and x, +1 are replaced by new randomly chosen x's. In a standard 
MF description, on the other hand, the notion of position is completely lost and one can, 
e.g., choose to replace the fitnesses of a certain number, K — 1, of other species selected at 
random, besides i m i n . This random neighbor (r.n.) model is the only one for which a MF 
treatment of avalanches could be set up so far p|JTT|. However, the lack of any meaning for 
distance in this MF is a quite strong limitation, to the extent that the very notion of SOC 
can be legitimately questioned |L2 . 



Avalanches correspond to sequences of mutations in which the minimal x species is always 
found among those resulting from genetic changes in previous stages, starting from a given 
ancestor mutation with x~ . = A. In the system the minimal x value does not exceed A 
for the whole duration of the avalanche. The probability P that an avalanche involves s 
mutations is expected to vary asymptotically as P(s) oc s~ T in the SOC state, in which for 
all avalanches A attains the value x c , the sharp threshold of the stationary x-distribution ||. 

So far, in models like the n.n. chain, r and similar exponents have always been discussed 
as bulk quantities |12 , ^3[ , i.e. considering statistics of avalanches starting everywhere within 



large, periodic systems. Compared to those in the bulk, a species at one end of an open chain 
(e.g., main predator, or basic level of microscopic life) has less species directly or indirectly 
connected to it. The paths through which dynamical correlations can propagate starting 
from an initial mutation on the boundary are also reduced. So, e.g., in semi-infinite geometry, 
boundary avalanches could be characterized by peculiar exponents, different from the bulk 
ones. Demonstrating boundary scaling in models like the BS one is a challenge, especially at 
the analytical level. Indeed, in the context of SOC with extremal dynamics, exact results are 
essentially limited to the above mentioned MF treatment [ |TT| , |T1| | . Consideration of boundary 
effects or other inhomogeneities clearly requires a meaningful notion of distance. We achieve 



this within a novel MF description of the BS model with n.n. interactions, generalizing the 
standard BP studied in probability theory M. 



The main scaling result for the random neighbor MF model is r = 3/2 PJII|. This r is 
consistent with MF BS avalanche dynamics being equivalent to a BP. An avalanche can be 
identified with a tree, where nodes represent species mutating within the avalanche. From 
each node as many branches depart, as there are species undergoing genetic change directly 
due to a mutation taking place at that node. The same species can act as node more than 
once within an avalanche. The complex structure of correlations of the BS model is simpli- 
fied in MF by assuming that, at each node, well defined, independent probabilities exist for 
all branchings compatible with the dynamical rules. Avalanches are generation trees, whose 
distribution in number of generating individuals, s, is given by P. The existing MF ap- 
proach clearly can not address exponents for diverging lengths, as defined, e.g., in a Landau 
approach to standard criticality. We introduce a characteristic length within MF through 
boundaries breaking translation invariance and leading to a position dependence of the BP 
description. Standard BP theory deals with the discrete transform P(z) = P{s)z s , on 
which the scaling of P(s) produces singular behaviour of the form 

P(z) = 1 -c(l -z) 1 " 1 +l.s.t. (1) 

for z — > 1~. In eq.(|l|) c is a suitable positive constant and the last term on the r.h.s. indicates 
regular or subleading singular terms. Without making reference to relative locations of 
the species along the chain, the standard BP assumes that well defined probabilities, Pi 
(i = 0, 1, 2, ...K), apply to the events in which a given species undergoing mutation triggers 
subsequent genetic changes, in the same avalanche, in % species, possibly including itself. 
Independence of branchings leads to validity of Watson's functional equation [|TJ 

P(z) = zG(P(z)) (2) 

with G(y) = p + Piy + P2IJ 2 + ••• +PkV K ■ Eq.(0) imposes a constraint on the p^s consistent 
with a singularity of the form (P. Such constraint reads G'(l) = 1 and automatically fixes 



c = J2/G"(l) and r = 3/2 as the only compatible exponent \TE\. This result for r is 



largely universal with respect to different choices of the parameters pi and only relies on the 
analyticity of G. A natural choice is Pi = x l c (l — x c ) K ~ l . In the r.n. model x c — 1/K 0, 
implying satisfaction of G'(l) = 1. Replacing x c by A < x c would amount to consider off- 
critical avalanches (G'(l) < 1), with x c — A playing the role of temperature-like field. Let us 
consider now a semi-infinite sequence of species on a chain. To each species is associated an 

integer coordinate j = 0, 1, 2, In a n.n. model the presence of the boundary requires to 

allow for a j-dependence of the avalanche size distribution: thus, Pj(s) or Pj(z) will describe 
avalanches starting at site j along the chain. This situation can still be analysed within what 
we call here inhomogeneous BP. Since, as a consequence of a mutation at j > 1, at most 
3 species can be further involved in the avalanche (K = 3 for the n.n. case), probabilities 
Po,Pi,P2, and ps will describe the possible outcomes of such a mutation. For convenience, 
and consistently with the above expressions of the p^s in terms of x c , one can further assume 
that with probability po, no further mutation takes place in the avalanche; with probabilities 
Pi/3 and P2/3 the avalanche propagates, respectively, in any one and any two of the species 
in the set {j — l,j,j + 1}; finally, p% is the probability that the avalanche involves all three 
species. In the MF spirit it is also sensible to assume j-independence for the p^s as long as 
j > 1. Of course, there should be different probabilities p\ for j = 0, where the boundary 
imposes p' 3 = 0. A possible choice made below is to assign p' = p + ^pi,p[ = \p\ + \pi and 
p' 2 = |p 2 +P3 at j = 0, again implying equivalence of j = and j = 1 with respect to single 
branch outcomes. 

With the above positions, eq.(2) is replaced by a full hierarchy of equations: 



p' 2 (P (z)P 1 (z) 



P j ( z ) = z (p + Pl(p j _ 1 ( z ) + P 3 ( z ) + 
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P j+ i(z)) + f(P j - l (z)P j+1 (z)+ 
P j (z)P j+1 (z) + P j - 1 (z)P j (z)) + 



P3 



(Pj-MPjWPj+tiz))) ;j > 1. (3) 



Pj(z) should converge to the bulk solution of eq.(^), for j approaching infinity. Thus, it 
is advantageous to adopt the following ansatz: 

Pj{z) = P{z) + A{z)e- qiz)j + l.s.t. (4) 

where q is an inverse length and P is the solution of eq.(^). As shown below, the assumed 
j-independence of A and q is consistent, as corrections to it would only involve subleading 
singular terms for z — > 1~. By substituting eq.(f|) into eqs.© one can deduce singular 
behaviours of Pq and q. For z — ► 1~, we expect A(z) ~ (1 — z) a and q(z) ~ (1 — z) 13 , with 
a and f3 suitable exponents. After substitution in eqs.(|]) for j > 1 one gets 

1 = |(l + 2coshg(^)) (G'(P(z)) + 

(5) 

Taking into account that P has the form ([j]) with r = 3/2, the leading singular terms in 
eq.(|) give: 

^- ) + l -A{z)=a{l-zfl 2 + l.s.t.- (6) 

where a = c of eq.(|I|). The same kind of substitution in the first of eqs.(S) leads to 

A(z) = a(l - z) 1 ' 2 - bA{z)q{z) + l.s.t. (7) 

with 6=1. Eq.(||) and eq.(|7|) determine both a and (3 above. In particular Pq(z) takes the 
form 

P (z) = P(z)+A(z) + l.s.t. = 1-6(1 -zf^ + l.s.t. (8) 

In general b = 2 (p' 1 /2+p'')-i anc ^ ^he resm ts @ and (H) make sense for J2w'i < 1- This 
condition is satisfied by our choice of p-'s, which further acquire the form p[ = ( K j rl ) x c(^ — 
x c ) K ~ 1 ~ i , if the pj's are expressed in terms of x c as discussed above. Thus, the threshold x c 



for the distribution of x values at the borders in the stationary state should be the same as 
in the bulk. According to eqs.([T|) and (g) 

Po(s) ~ s~ 7/4 (9) 

Thus, in our MF description the BS SOC state is characterized by a boundary scaling with 
an exponent t' = 7/4 different from the bulk one. Boundary avalanches of course suffer 
more rapid extinction and their distribution decreases faster for large s. It is interesting 
to note that, by exploiting analogies with magnetic systems, r' = 7/4 has been predicted 
recently within a MF approach to border avalanches in Abelian Sandpile Models (ASM) 



with Dirichelet boundary conditions [fLq| . This lends further support to the idea that in 
ASM a BP description underlies the statistics of avalanches in the MF limit, for which also 
t = 3/2 is expected [|17]]. By a numerical approach one can also identify r' ~ 7/4 for MF 
avalanches of the earthquake model of ref. || ||18|| , confirming an underlying BP also in this 
case. A further consequence of eqs.(|6|), (0) is the singularity: 

q(z) ~ (1 - zff\ (10) 

Thus, the penetration length of the border disturbance, q~ l , diverges for z — > 1~. In 
MF treatments of inhomogeneous equilibrium models, quantities like q~ x show the same 
divergence with temperature as classical correlation lengths. By interpreting z as a standard 
fugacity for a polymer, one deduces from eq. (|I~0"D a correlation length exponent v = 1/4. This 
is indeed the classical v of branched polymers |19[. Of course the definition of v for a SOC 
system requires one to identify physically meaningful parameters describing the approach or 
the departure from criticality. For BS avalanches such a parameter is the temperature-like 
deviation x c — A. By introducing A-dependent p^s and p-'s in our equations, the result fllUD 
can be converted into g(A, z — 1) ~ (x c — A) 1 / 2 , which implies v — 1/2. Remarkably enough, 
this is the classical v exponent expected for ASM [IT]. This and the above mentioned 
coincidence of r' strongly support the idea that BP fully underlie also the MF description 
of ASM avalanches. 



In order to identify boundary scaling beyond MF, we performed systematic simulations 
with open, n.n. BS chains of different lengths (N < 10 3 ). First we verified that the 
distribution of boundary x's in the stationary state is essentially unaltered with respect to 
that of the periodic, bulk case, and displays the same sharp threshold at x c = 0.665 ± 0.015 
0. This coincidence is fully consistent with our choices of the p'^s in the MF approach. By 
selectively sampling avalanches starting near the boundaries or in the interior of the chains, 
we extrapolate r' = 1.25±0.01 (see FIG. 1). This value is clearly different from the bulk one 
t ~ 1.08 fl2|]. So, also in the n.n. model boundary avalanches have a probability decaying 



more rapidly at large s, than in bulk. A further characterization of boundary scaling is given 
by the distribution of first return times of activity (x taking the minimum value) at the same 
boundary site. These times are distributed as t~ T f irst , with 7f irst = 1.35 ± 0.01, different 
from the bulk value Tf irst = 1.58 |13). By recording the times of all subsequent returns of 
activity one can also obtain a distribution oc t~ r «», with T' all = 0.65 ± 0.01, again distinct 
from T a u — 0.42 in bulk fI5| . Such boundary exponents are consistent with a scaling relation 



T 'first + T 'aii = 2, already satisfied in the bulk [[□[]. Since the validity of such relation should 
not depend on the position considered along the chain, the above consistency is further 
indication of the good quality of our determinations. Data concerning these exponents are 
shown in FIG. 2. 

We conclude that at the boundaries activity has a different pattern of intermittency. 
First returns are shifted towards longer time scales. On the other hand, once the boundary 
has been reached, activity remains more easily trapped there, giving rise to concentrated 
sequences of returns. In applications of the BS model, the choice of a more or less reg- 
ular network of interactions remains to some extent arbitrary, and should not matter for 
universal properties, unless the long-range limit of a r.n. model is assumed. However, the 
distinction elucidated above between bulk and boundary species, appears to have important 
consequences, affecting the universal scaling features of evolution. Thus, boundary scaling 
offers additional, deeper insight into the properties of biological models and widens the con- 
text of their possible comparison with paleontological data. Summarizing, we showed here 
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that within the framework of punctuated equilibrium there exists a well defined boundary 
scaling in addition to the bulk one. At the MF level this scaling can be analysed exactly 
within a generalization of BP theory, which considerably extends previous classical descrip- 
tions of the BS model, and directly focuses on its relation with other models. In particular 
also ASM with Dirichlet boundary conditions fall fully in the MF universality class of our 
BP. Also in the n.n. case our results show the existence of new scalings which make the 
notion of species at the ends of a chain meaningful in a universal sense. 



9 



REFERENCES 

[1] T. Harris Theory of Branching Processes, Springer, Berlin (1963); W. Feller An Intro- 
duction to Probability Theory and its Applications, Wiley & Sons, New York , (1971). 

[2] P. Bak, C. Tang and K. Weisenfeld, Phys. Rev. Lett. 59, 381, (1987). 

[3] P. Alstr0m Phys. Rev. A 38, 4905 (1988); S. Zapperi, K. B. Lauritsen, H. E. Stanley 
Phys. Rev. Lett. 75, 4071 (1995) 

[4] L. Pietronero, A. Vespignani and S. Zapperi Phys. Rev. Lett. 72, 1690 (1994). 

[5] Z. Olami, H. J. S. Feder and K. Christensen Phys. Rev. Lett. 68, 1244 (1992); S. Lise, 
H. J. Jensen preprint. 

[6] A. E. Scheidegger Bull. Assoc. Sci. Hydrol. 12 15 (1967); H. Takayasu J. Stat. Phys.bf 
65, 725 (1991) 

[7] P. Bak and K. Sneppen Phys. Rev. Lett. 71, 4083 (1993). 
[8] H. Flyvbjerg, K. Sneppen and P. Bak, Phys. Rev. Lett. 71, 4087 (1993). 
[9] N. Eldredge, S. Gould, Nature (London) 332, 211 (1988). 
[10] D. M. Raup Science 231, 1528 (1986). 

[11] J. de Boer, B. Derrida, H. Flyvberg, A. D. Jackson, T. Wettig Phys. Rev. Lett. 73, 906 
(1994). 

[12] J. de Boer, A. D. Jackson and T. Wettig Phys. Rev. E 51, 1059 (1995). 

[13] S. Maslow, M. Paczuski, P. Bak Phys. Rev. Lett. 73, 2162 (1994). 

[14] As an exception, see S. Boettcher and M. Paczuski Phys. Rev. Lett. 76, 348 (1996) 

[15] A condition analogous to our G"(l) = J2i Wi — 1> is also at the basis of the RG approach 
to sandpiles of ref. [4], where it enforces energy balance. 



10 



[16] A. L. Stella, C. Tebaldi and G. Caldarelli, Phys. Rev. E 52, 72 (1995). 
[17] C. Tang and P. Bak Journ. Stat. Phys. 51, 797 (1988). 
[18] S. Lise, A.L. Stella to be published. 

[19] A. B. Harris, T. C. Lubensky Phys. Rev. B 24, 2656 (1981). 



11 
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FIG. 1. Q(s) = f* max P(s')ds' is the integrated distribution; the fitting form is As 1 ""' + C with 
t' = 1.25 ±0.01 (N = 1000). 
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FIG. 2. first and all return times probabilities at boundary sites. Statistics refer to 10 9 muta- 
tions in the whole chain. 
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